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This research paper investigates current experimental designs and regression models for calibrating 
internal wind tunnel force balances of non-monolithic design. Such calibration methods are necessary for 
this class of balance because it has an electrical response that is dependent upon the sign of the applied forces 
and moments. This dependency gives rise to discontinuities in the response surfaces that are not easily 
modeled using traditional response surface methodologies. An analysis of current recommended calibration 
models is shown to lead to correlated response model terms. Alternative modeling methods are explored 
which feature orthogonal or near-orthogonal terms. 

Nomenclature 

a t = intercept 

Fj = factor setting 

R t = model response for a particular degree of freedom 

Txy = Pearson’s correlation coefficient 

x = independent variable 

Z = indicator variable 

BBD = Box Benken Design 

CCD = central composite design 

LaRC = Langley Research Center 

LSR = least squares regression 

RSM = response surface methodologies 

SYS = single vector system 
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TASK= type of wind tunnel force balance manufactured by the TASK corporation 


VIF = variance inflation factor 


I. Introduction 

A typical wind tunnel internal balance is a six degree of freedom force and moment transducer, capable of 
measuring an aerodynamic normal force, side force, axial force, pitching moment, yawing moment and rolling 
moment by monitoring deflection of structures and their resultant strain. The balance is said to be “internal” 
because the forces and moments are measured within the wind tunnel model. 



Figure 1 Typical Six Degree of Freedom Internal Wind Tunnel Force Balance 

Figure 1 is a typical six degree of freedom internal wind tunnel force balance used for measuring aerodynamic flight 
characteristics. The metric (left) end of the balance connects to the aircraft model, while the non-metric (right) end 
is where the balance is constrained to the sting. There are three measuring sections within such a balance; two cage 
sections, and an axial section. The cage sections are used to measure all forces and moments except axial force. 

The axial section measures one force, axial force. Typically these sections are machined out of one solid block of 
material (called a monolithic design) to avoid problems with mechanical hysteresis. 

The type of balance under consideration for this research has two measuring sections, the cage section and 
the axial section, each machined individually and fastened together mechanically; hence the balance is a non- 
monolithic structure. Mechanical joints exhibit a different response behavior than that of a monolithic balance. 
Monolithic balances display a continuous response output when the applied load transitions from positive to 
negative or from negative to positive in sign. In contrast, the output of a non-monolithic balance may show a spike 
at the zero load point and may also have a different slope on either side of the zero load point as shown in figure 2. 
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Figure 2 Non- Monolithic Force Balance Response due to Applied Force 

The response behavior across either side of the zero load point may not be limited to just a change in slope; 
it could also affect a change in quadratic curvature or other parametric curvature. The change in response across 
the zero load point may also be different for each sub-domain of the entire design space. This means that a response 
may be dependent upon the polarity of all forces and moments. The AIAA recommend practices document (AIAA 
R-09 1-2003) states, “. . .it is not uncommon for the load/output relationship of balances, especially those of multi- 
piece design, to exhibit some dependency on the sign of the strain in the measuring element. This asymmetry results 
in the need to determine and use different calibration coefficients according to the sign of the force or moment 
acting on the bridge, in order to achieve the best accuracy from the balance [1 ].” The “asymmetry” that is referred 
to is a challenging phenomenon to model because it cannot be captured using traditional balance calibration methods 
and its behavior is not fully understood. In the following paragraphs, it is shown that the current method suggested 
in R-09 1-2003 does not sufficiently correct for this problem. In order to improve the accuracy of non-monolithic 
balances this modeling challenge must be further investigated. 

IT. Current Practice 

One of the goals of the recommend practices document (AIAA R-09 1-2003) is to formalize the calibration 
process, leading to a universal calibration matrix form that can be used for comparison between calibrations 
performed at different institutions [1]. The recommended practice for calibrating non-monolithic balance suggests 
building a single regression model for each of the six degrees of freedom, as opposed to building multiple models in 
individual sub-domains for each degree of freedom. The former is preferred over the latter since it provides a 
simple, compact method for exchanging calibration matrices within the engineering community. The recommended 
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practice suggests using an extension of the traditional response surface model with the inclusion of absolute value 
terms in order to model asymmetric behavior. The 97 term model is shown below in equation 1. While this is the 
recommended standard model, it has not been fully implemented in all calibration institutions. 


7=1 
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In the suggested model, terms that include an absolute value will act as a modifier to the similar effects that 
lack the absolute value. The model attempts to create terms that will adjust the response surface depending upon the 
polarity of the applied forces and moments. Regardless of whether or not this modeling strategy properly captures 
asymmetric behaviors, the 97 term model suffers from poor model quality metrics that are inherent to the model 
formulation. These metrics quantify the collinearity between model terms and the prediction variance of future 
observations, both of which suffer from over-parameterization. The recommended practice suggests discarding 
terms in the model when necessary; however there are certain terms in this model that should not coexist. To 
illustrate for example, a simple comparison is made between a quadratic response and the absolute value of a cubic 
response, shown in figure 3. By visual inspection of figure 3 it is obvious that the quadratic function and absolute 
value of a cubic function have similar modeling effects. 
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Figure 3 Comparison between absolute value of cubic and quadratic terms 
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Another similar example involves the quadratic term and the absolute value term. In figure 4 it appears the 
degree of collinearity between terms is less than that of figure 3. The deviation of the quadratic curve from the 
absolute value response may be great enough to lead one to believe that these terms are capturing significantly 
different behaviors, thus, justifying the inclusion of both in the overall model. To determine the degree of their 
ability to model unique behavior, or worded differently, the severity of their collinearity, model quality metrics must 
be assessed. 

III. Model Quality Metrics 

The "design" is the integral part of a designed experiment, from which most of the statistical quality metrics 
of the experiment are derived. The phrase "design," also referred to as "design matrix" was coined in the fields of 
Design of Experiments (DOE) and Response Surface Methodology (RSM) to refer to the combinations of settings 
(levels) of all independent variables. A model matrix contains information about the test schedule for the individual 
calibration points as well as information about the model form. For demonstration, a calibration design matrix, 
model matrix and regression model are shown in equation (2.1), (2.2) and (2.3) respectively. 

NF SF 
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Each row of the calibration design, D, represents a different run in the calibration experiment. The number of 
columns in a model matrix depends upon the number of coefficients in the regression model. The predicted bridge 
output, y, models bridge output as a function of applied forces. 

A. Pearson Correlation Coefficient 


For a given model matrix, the Pearson correlation coefficient r xy describes the degree of the collinearity 
between columns in the model matrix [2]. 


r xy ~~ 
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It can be shown that the range of r xy is: 


-1 <r X y< 1 


Generally values of r xy near 1 indicate a strong positive linear association between x and y whereas values of r xy 
near -1 indicate a strong negative linear association. Values of r xy near 0 indicate little or no linear association 
between x and y . Computing the Pearson correlation coefficient between each column in a model matrix leads to a 
symmetrical correlation matrix R shown below. 
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B. Variance Inflation Factors 


Variance inflation factors (VIF) are commonly used in place of correlation coefficients because of their 
more compact description of multicollinearity in a design [3]. The VIF measures how much the variance of the 
model is increased by the lack of orthogonality in the design. True orthogonality guarantees that each regression 
coefficient estimate is independent of any other. Each term in the model (excluding the intercept) has a calculated 
VIF, where a VIF of one indicates orthogonality. As a VIF increases, the collinearity of that term with others in the 
model also increases. The VIF’s of a design are computed as the diagonal component of the inverse of the Pearson 
correlation matrix. A general rule of thumb is that a VIF greater than ten is an indication that multicollinearity may 
be unduly influencing the least squares estimates. Some believe this rule may be too lenient and suggest that a VIF 
should not exceed five. A limitation to the use of VIF’s is that they cannot distinguish between several simultaneous 
multicollinearities. When high VIFs are encountered it is useful to refer back to the Pearson correlation matrix to 
determine which factors in particular are correlated [3]. 


IV. Case Study of a Non-Monolithic Balance Calibration Design 

The load schedule of a non-monolithic balance calibration, carried out at NASA LaRC has been provided. 

Using this data, the model matrix for the 97 term model was reconstructed which was used to compute the 
aforementioned model quality metrics. This sample load schedule varied two factors at a time. One factor was held 
constant while another swept its respective domain. This procedure was repeated until all possible combinations 
were executed. The load schedule consisted of 1,316 runs, resulting in a model matrix consisting of 1,316 rows and 
97 columns. From the model matrix, the following variance inflation factors were computed in table 1. 


Table 1 Mean VIF for individual model terms 


Term 

Approximate Mean VIF 

H=iblijFj 

132.9 

Yj=i 

148.0 

2?=iCl ijFj 2 

1,241.0 


1,219.5 

£? =1 S?=,- +1 C3 iJik FjF k 

1.3 

Yj=i Hk= j+1 \PjFk I 

2.9 

Z%iIZ=i + ic5ij*Fj\F k \ 

1.9 

Ytj=\Ytk=j+l c ^i,j,k \Fj\F* 

2.0 

Xj=i dlij F? 

1,084.2 

X%idluW\ 

1,086.2 
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By inspection of the Pearson correlation matrix the following terms are most highly correlated (shown in table 2). 


Table 2 Mean Pearson correlation coefficients for the most highly correlated terms 


Term* 

Term y 

Approximate Mean R xy 

£y=i b\ijFj 

X7=ic2ijFj\Fj\ 

0.97 

E?=i c Vv 2 

Yj=ib2ij\Fj\ 

0.96 

S?=iW ijFj 

Y.%td\ Uj Ff 

0.91 

Ey=i dlij |F/| 

Y!j=ib2ij\Fj\ 

0.91 

wi 

23 Li cXi'jFf 

0.99 

Z?.i ilyf? 

27=1 ^i.jpj\Fj\ 

0.99 


The consequence of highly correlated model terms is large prediction variance estimates in the design 
space. The variance of a predicted value is a function of the model matrix and is calculated as follows [3,5]: 

Var\y(x )] = x^ f (X f X)~ 1 x^ ■ a 2 (33) 

where x^ is a function of location in the design space and also a function of the model, X is the model matrix, and 
<j 2 is the mean squared error. For comparative purposes, it is commonly assumed that a 2 equals one, therefore in 
terms of standard error of prediction equation 33 becomes [3,5]: 

JVar\y(x)\ = (34) 

The multicollinearity diagnostics presented in table 1 and 2 are derived from the model matrix, which 
depends on the experimental design and model formulation; therefore, the cause of the large variance inflation 
factors can be attributed to both the experimental design and model form. Decoupling the amount of collinearity 
attributed by each is difficult because no experimental design exists for the 97 term model that yields orthogonality. 
Therefore, it is challenging to state that a certain amount of collinearity is due to the experimental design, because 
there is no benchmark to compare it to. Later in the paper, it is found that multicollinearity in a newly proposed 
model is completely attributable to the 1,316 run experimental design and not to the model form. This analysis 
highlights a fundamental problem with the 1,316 run design and demonstrates the need for careful consideration of 
the experimental design when employing different models. Regardless of the source multicollinearity, the variance 
inflation factors presented in table 1 are unacceptably large. This reflects a problem with both the 1,316 run 
experimental design and the recommend practices 97 term model. 
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V. Alternative Models 


Four alternative methods for modeling the response of non-monolithic force balances are presented here. 
The first approach suggests building individual models in separate sub-regions. The second is to use a pure cubic 
model to capture the entire response with a global model (as per the standard). The third is an absolute value model 
approach, while the fourth method uses an indicator variable scheme to fit a global regression model. 

Both the independent model and indicator variable approach suggest partitioning the entire design space 
into sub-spaces. Due to typical force ratios inherent to balance design, normal force is the highest load exerted on 
the joints of the multi-piece balance and thereby produces the most prevalent asymmetric behavior. For the 
independent models and indicator variable methods, the design space will be partitioned into sub-spaces dependent 
on the direction of the applied normal force. Regardless of how the other five components are loaded, when normal 
force is negative, this will correspond to a sub-region, and when normal force is positive this corresponds to another 
sub-region. 


A. Independent Models Approach 

The first approach suggests using separate independent regression models and experimental designs for 
each sub-region. Designs exist such that the model shown in equation 35 can be fit while having all variance 
inflation factors approximately equal to one. A recommended design to implement in each sub-space is a central 
composite design (CCD) which is commonly used to build the second order RSM model as shown below [3]: 

n n n n 

bl i,j F j + ^ b 2 i,j 

7=1 7=1 7=1 fc=7+l 

By building separate models, the quality metrics for each model can be calculated individually. If a CCD is chosen 
for each subspace, it can be shown that the quality metrics of each CCD model will be ideal. That is, all variance 
inflation factors are one, and there is no collinearity present in the model [3]. One disadvantage of this approach is 
that it requires some form of a conditional statement in the data acquisition software so that the proper regression 
models are used depending on the applied forces. Another disadvantage is that it has fewer degrees of freedom for 
the lack of fit analysis as compared to the indicator variable approach presented later. Since the independent 
models approach would perform a least squares regression (LSR) in one sub-region, only the design points used in 
that sub-region could be utilized in the lack of fit analysis. In contrast, the indicator variable approach would 
perform LSR over multiple sub-regions, thus providing more degrees of freedom for estimating lack of fit. 


*♦1 z b3 iiJik FjF k (35) 
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Since the late 90’s the RSM model (equation 35) fit with a 2 6_1 spherical CCD using the Single Vector 
System (SVS) has become the standard for calibration practice at NASA LaRC [1,6]. Adapting this method to a 
non-monolithic balance may not be the most innovative approach but it is arguably one of the most reliable. 

B. Pure Cubic Model Approach 

The pure cubic model attempts to capture the asymmetric response behavior with a global model over the 
entire region of operability. Previous research has been conducted by Draper that has led to variance optimal 
rotatable designs for the pure cubic model[7,8]. The design has three levels made up of a central composite and a 
Box Benken Design. The first orbit contains factorial points while the second orbit has the axial points of the CCD 
and design points of the BBD. A third and final orbit makes up the outer most axial points [9]. While Draper’s 
design uses axial points necessary for a rotatable design, a spherical design may be more practical when extending 
Draper’s design into six factors. It can be shown that Draper’s design in six factors (228 runs) can support the pure 
cubic model shown in equation 36 with variance inflation factors below three. 

n n n n n n n 

Ri = &i + ^ 

7=1 7=1 7=1 k=j+l m j=m+lk=m+j+l 

n 

+ ^Z>5 y F/ (36) 

7=1 

The pure cubic model of equation 36 was not developed for the specific application of classifying the asymmetric 
behavior of a non-monolithic force balance. Instead, it is proposed as a heavily parameterized model that does not 
suffer from overwhelming multicollinearity. This model takes the approach of the 97 term model by including as 
many parameters as possible in order to provide the greatest flexibility of shape for the response surface estimation. 
The use of equation 36 allows for use of standard regression model reduction methods without the ambiguity 
associated with the 97 term model. For these reasons, the pure cubic model is a good candidate for modeling multi- 
piece balance behavior. 

To demonstrate the deficiency of the 1,316 run experimental design presented in the case study, VIF’s of 
the pure cubic model are compared using the 1,316 run experimental design and Draper’s 228 run experimental 
design, as shown in table 3. 


^ ' blj'jFj F ^ ^ ' b3ij ik FjF k + ^ ^ ^ ' b^i.mJ.kRrnRjRk 
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Table 3 Comparison of Mean VIF’s for the Pure Cubic Model 


Term 

1,316 Run Schedule 

228 Run Schedule 

27-1 M u f j 

4,961.3 

2.0 

17-1 MyF/ 

1.1 

1.0 

Z7=iIZ=j + ib3 iJJC FjF k 

324,389.0 

1.0 

lln l]=m + l IZ=m +J+ l bhm,j, k F m FjF k 

241,296.3 

1.0 

bS u F? 

6.8 

2.0 


From table 3 it is clear that Draper’s experimental design demonstrates superior orthogonality between model terms. 
Additionally, Draper’s design benefits from a time and cost perspective by using only a fifth of the number of runs. 
Since orthogonality is preserved using one experimental design and is lost using another, the conclusion can be 
made in this case that the experimental design is the root cause of the multicollinearity and not the model form. 

C. Absolute Value Model Approach 

The third approach fits a global model using similar parameters to that of the 97 term model. This 
approach models main effects, absolute value of the main effects and two factor interactions. The model, referred to 
as “the absolute value model” is shown below: 

n n n n n 

Ri = a i + '^ bl u Fj + b2 u | Fj | + £ b3 u Fj\Fj \ + £ £ b4 iiJik FjF k (37) 

7=1 7=1 7=1 7=1 k=j+l 

It is a reduced model compared to the 97 term model, having only 34 model terms. The strengths of this model lie 
in its ability to capture first order terms in the presence of hysteresis while maintaining the orthogonality of all 
terms. To illustrate the capabilities of the absolute value term model, consider a purely linear single factor 
response with a discontinuous slope that changes at the origin, shown in figure 7. 



Figure 7 Asymmetric Slope Modeling Comparison 
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The limitations of the quadratic model become clear; the model is fit so that the quadratic curvature accounts for the 
slope change. In balance calibration, first order effects are usually the most statistically significant effects. Hence, 
using an uncorrelated, absolute value model may be the best approach to capture this dominant effect particularly if 
there is no requirement for correcting the response for higher order interactions. 

To demonstrate another deficiency of the 1,316 run experimental design, VIF’s of the absolute value model 
are compared using the 1,316 run experimental design and Draper’s 228 run experimental design , as shown in table 
4. 


Table 4 Comparison of Mean VIF’s for the Pure Cubic Model 


Term 

1,316 Run Schedule 

228 Run Schedule 

SjLi blijFj 

16.7 

4.0 

ZUbliM 

1.3 

1.2 

lU b hjFj\Fj\ 

1.2 

4.0 

S?=i512=, + 1 b4 Uk F,F k 

15.9 

1.0 


D. Indicator Variable Approach 

Consider now a method to model sub-regions separating asymmetric behavior using a global model with 
indicator variables [4]. Indicator variables are commonly used to model different processes that are independent of 
each other. Classically, the independent processes span the same domain where the indicator variable has no spatial 
dependence. In this approach, the indicator variable is used as a spatially dependent categorical variable that 
changes depending on the sign of the applied normal force. Similar to the independent models approach, there is a 
design in each sub-region. In this case, the LSR is performed over the entire design space so that there are twice as 
many design points as the independent models approach but typically fewer than twice as many parameters. Since 
the increase in design points is greater than the increase in model parameters, extra degrees of freedom are awarded 
to the indicator variable approach. Additionally, this approach gains preference over the independent models 
approach since it can be presented in a single compact model that can be easily compared across institutions. 
Furthermore, the examples provided for this method display zero collinearity. The first example provided is a 
simple case to demonstrate the functionality of the indicator variable method. The single factor model shown below 
contains an intercept, (3 0 , a main effect, Pix b and an adjusting term to the main effect, x ¥ 1 . 

Y = p 0 +/3 1 X 1 + 'V 1 Z 1 X 1 (38) 
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When using indicator variables, a guideline for choosing the proper indicator setup is that there should be one less 
indicator variable than there are groups or sub-spaces being modeled. For this case there are two sub-regions, 
therefore one indicator variable, Z, is used [4]. The indicator variable is dependent upon the applied normal force, 
X l9 and is setup as follows: 

Z = -1, ifX 1 < 0 (sub-region# 1) (39) 

Z = 1, if Zi > 0 (sub-region#2) (40) 

The coefficient estimates, jS ± or W l9 alone, cannot model the response of either region, but when is subtracted 
from p 1 this corresponds to the response of sub-region number one, and when X ¥ 1 is added to p 1 this corresponds to 
the response of sub-region number two. For example, consider the following single factor response: 

? = 0.25^-0.15 Z t X ± (41) 

Where if Z = — 1 (sub-region# 1) the model becomes: 

? = 0AX 1 (42) 

If Z = +1 (sub-region#2) the model becomes: 

? = 0.1*! (43) 

Next it is shown that the previous example can be expanded into higher dimensions while maintaining 
orthogonality. Since it is assumed only the slope of p 1 changes with the direction of the applied X l9 fi 2 in this two 
factor example remains constant across sub-regions. The former model is simply extended to incorporate a second 
factor, X 2 , as follows: 

Y = P o + p 1 X 1 + fi 2 X 2 + ^ZA (44) 

A representative response surface modeling by the equation above is shown below in figure 8. 
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Figure 8 Response Surface with Asymmetric Slope in Single Direction 

When choosing a design for this type of model, it seems appropriate to build a traditional factorial design in 
each sub-region. It should also be noted that the Z column is omitted from the design matrix when performing LSR 
or calculating quality metrics such as the correlation matrix or variance inflation factors, since the intercept is 
constant across the entire domain. The resulting correlation matrix is a three by three identity matrix, while the 
VIFs are a row vector of ones, implying perfect orthogonality. 

Using the same procedure it is possible to extend the model to include six factors, with one linear adjusting 
term shown below in equation 45. 

6 

P = 4^1*1 (45) 

1 = 1 

The next example extends equation 45 further by introducing a two factor interaction term shown below: 

? = 00 + 01*1 + 02*2 + ^ 1 * 1*1 + 012 * 1*2 (46) 

Similar to the assumption made about the slope, /? 2 , it is also assumed that /? 12 is constant across the entire design 
space. A response surface representative of this model is shown in figure 9. 


14 



Figure 9 Asymmetric Slope in Single Direction with Two Factor Interaction Effects 

The same design is used for this model as was used for the previous example, with the addition of the two factor 
interaction column, /? 12 . Once again, orthogonality is preserved resulting in VIF’s of one for all model terms. 

Limitations to the indicator variable approach arise when trying to model higher order terms such as a 
quadratic along with the adjusting term to the main effect. The adjustment term is analogous to the absolute value 
term in that they both effectively model the same phenomena and are correlated. Previously in this paper, it was 
shown analytically that quadratic and absolute value terms could not be orthogonal in a design. Similarly, it can be 
shown that a quadratic term and the adjustment to the main effect term cannot coexist orthogonally. 

VI. Conclusion 

Either one of the four newly proposed modeling approaches is a viable alternative to the recommended 
practices 97 term model. Comparisons made between all approaches, including the 97 term model, used metrics 
that are calculated in the design planning and model formulation stage. Significant insight into model quality can be 
obtained in the planning phase of a calibration experiment. Careful consideration should be placed on coordinating 
the experimental design with the model formulation. 
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